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Abstract 

We present a set of efficient techniques in first-principles electronic-structure calculations uti- 
lizing the real-space finite-difference method. These techniques greatly reduce the overhead for 
performing integrals that involve norm-conserving pseudopotentials, solving Poisson equations, 
and treating models which have specific periodicities, while keeping a high degree of accuracy. 
Since real-space methods are inherently local, they have a lot of advantages in applicability and 
flexibility compared with the conventional plane-wave approach, and promise to be well suited for 
large and accurate ab initio calculations. In order to demonstrate the potential power of these 
techniques, we present several applications for electronic structure calculations of atoms, molecules 
and a helical nanotube. 

PACS numbers: 31.15.Fx, 02.70.Bf, 71.15.-m, 73.22.-f 



*Electronic address: ono@upst.eng.osaka-u.ac.jp 



1 



I. INTRODUCTION 



So far, a number of methods for first-principles electronic-structure calculations imple- 
mented entirely in real space have been proposed . These have 
desirable properties compared with the usual plane- wave and/or the linear combination of 
atomic orbital approaches: (i) Since all of the calculations are carried out in real space, it 
is easy to incorporate localized Wannier-type orbitals, which are localized in a finite region 
required for the realization of O(N) calculations, into the algorithm, (ii) A technique utiliz- 



ing a real-space double-grid 11 1 is available within the real-space finite-difference (RSFD) 
formalism, where many grid points are put in the vicinity of nuclei, so that the integrals 
involving rapidly varying pseudopotentials inside the core regions of atoms can be calcu- 
lated with a high degree of accuracy, (iii) In order to improve the calculation accuracy, 
the grid spacing should be narrowed, the procedure of which is simple and definite. Even 
more important is that (iv) boundary conditions are not constrained to beperiodic, e.g., a 
combination of periodic and nonperiodic boundary conditions for surfaces 12, ^ and wires 
and twist boundary conditions for helical nanotubes are applicable. This aspect (iv) 
is significantly advantageous to electron-conduction calculations, because an nonperiodic 
boundary is indispensable for the direction in which electrons flow 3]. 

In this paper, we present a set of efficient techniques for the RSFD approach. In order to 
demonstrate the applicability of these techniques, we calculate the electronic structures of 
atoms and molecules by the RSFD approach using isolated boundary condition and that of a 
(5, 5) carbon nanotube utilizing the twist boundary condition. In particular, the real-space 
calculations with the twist boundary conditions are appropriate for helical nanotubes, since 
the translational repeat distance leads to extremely large supercells, although the nanotubes 
can be constructed with one- dimensional translational symmetry. The rest of this paper is 
organized as follows: in Sec. II, III, and IV, we give in detail the efficient procedures for real- 
space calculations together with examples, and in Sec. V, we conclude with a discussion 
on the future direction of the RSFD calculations. The discretization of the Kohn-Sham 
equation in the RSFD formalism is addressed in Appendix. 



2 



II. TIMESAVING DOUBLE-GRID TECHNIQUE 



A. Essential feature of the timesaving double-grid technique 

In the RSFD formalism, wave functions, electronic charge density, and potentials are all 
represented on the discrete grid. In general, the discretizations of them are not invariant 
under uniform translations of the system respect to the position of the grid, which leads 
a serious problem in practical simulations: the total energy varies unphysically depending 
on the relative positions of grid points and the nucleus. This problem can be avoided by 
reducing the grid spacing; however, a fine grid may require so many points as to result in 
a substantial increase in computational effort. With this as a background, the timesaving 
double-grid technique was introduced Q to circumvent this problem easily. 

We here extend this technique so that it can be applied to local parts of pseudopoten- 
tials as well as nonlocal pseudopotentials. The double-grid consists of two sorts of uniform 
and equi-interval grid points, i.e., coarse and dense ones, depicted in Fig. 1 by "x" and 
respectively. The dense-grid region enclosed by the circle is the core region of an 
atom that is taken to be sufficiently large to contain the cutoff region of pseudopotentials. 
Throughout this paper we postulate that the wave functions and the charge densities are 
defined and updated only on coarse-grid points, while pseudopotentials are strictly given on 
all dense-grid points in an analytically or numerically exact manner. In the case of using 
the pseudopotentia! parameterized h y Bach** et al Q and Trouiher and Martins Q 
with the Kleinman-Bylander nonlocal form [17|], the ionic pseudopotential term in the total 
energy is defined as 



/ ,C A '(r - R')p(r)dr + [ v%"(r - R-)p(r)dr + & - ,9l ~' 



where M is the number of occupied states and 

9L,i= I vL(r - R s )Ur)dr. (2) 
Jn 

Here, p(r) is the electron density, and R s , v^J d ' s (r), and v s l °j. t,s {r) represent the position of 
the s-th nucleus, the local pseudopotential which varies rapidly around the nucleus (hard lo- 
cal part), and the local pseudopotential which changes gently (soft local part), respectively. 



The procedure of decomposing the local pseudopotential is explained below. In addition, 
v im( r ) = vf ( 7 *)V'to S ( 7 ")) where vf(r) and ipf^ s (r) are the nonlocal parts of the pseudopoten- 
tials and the pseudowave functions used to prepare the pseudopotentials, respectively. The 
numerical integrations related to v^ d)S {r) and vf m (r), which vary sharply in the vicinity of 
nuclei, are performed on the dense grids using the timesaving double-grid technique, whereas 
the integration related to the gently behaving part v^^ir) is implemented on the coarse 
grids. It is noteworthy that v^ d,s {r) and vf m (r) must vanish outside the core region of 
the pseudopotential; otherwise the calculations of Pulay forces, which are computationally 
demanding, are required in molecular-dynamics simulations. While the value of the nonlo- 
cal part is zero outside the core region, the value of the local part, which is a long-range 
Coulomb potential, is not zero. Thus, in the case of using the pseudopotential parameterized 
by Bachelet et al. [3], v^ d,s (r) is so set by the following procedure as to be zero outside 
the core region: 



v Har d ,s {r) = _ zs , ^vyp^ _ I , (3) 
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where a\ > a| and erf(x) is the error function (or probability integral). 

Let us first consider the inner products between wave functions ip and nonlocal parts of 
pseudopotentials vf m (see Fig. 2). Here, we assume a one-dimensional case for simplicity. 
The values of wave functions on coarse grid points (o) are stored in a computer, and the 
values on dense-grid points (•) are evaluated by interpolating them. The known values 
of pseudopotentials both on coarse- and dense-grid points (o) are also shown schematically. 
Then, from Fig. 2(a) one can see that only the values on coarse-grid points are so inadequate 
that the inner products can not be accurately calculated; the errors are mainly due to the 
rapidly varying behavior of pseudopotentials. On the other hand, Fig. 2(b) indicates that the 
inner products can be evaluated with a high accuracy, if the number of dense-grid points is 
taken to be sufficiently large and also if the wave functions on dense-grid points are properly 
interpolated from those on coarse-grid points |l8{ |. 

There are several interpolation methods for wave functions, among which the simplest is 
linear interpolation. In this case, the wave functions ipj = ip(xj) on dense-grid points Xj are 
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/d/2+R s nN core 
v! m {x)il>{x)dx « V v s lm ^jh dens , (6) 
-<2/2+R s „•_ „at 



interpolated from tyj = ip(Xj) on coarse-grid points Xj as 

^ = ^ + (5) 

where h is the grid spacing of the coarse-grid points. The inner product is assumed to be 
accurately approximated by the discrete sum over the dense-grid points, i.e., 

-d/2+R s nN CO r e 
-d/2+R<> , _ y 

J — lllvcore 

where t>f ■ = vf m (xj — R s ), d is the "diameter" of the core region, 2N core + 1 (2nN core + 1) 
is the number of coarse (dense)-grid points in the core region, hd ens is the grid spacing of the 
dense-grid points, and n = h/hd ens , i.e., n — 1 is the number of dense-grid points existing 
between adjacent coarse-grid points. Now, substituting Eq. (JSJ) into the right-hand side of 
Eq. (JOJ), we have 

/d/2 + R 3 N core 
v s lm {x)^{x)dx « w L,J^jh, (7) 

■d/2+R" r_ M 



where 

a \^ h — \x n j + k — Xj\ s 

W lm,J — 2-^1 nfo V lm,nJ+k- \°) 

k=—n 

Similarly, the inner product of the hard local part and electron density is given by 

/d/2+R s nN core 
hard,s / \ / \ j ^, \ ^ hard.s r 
V loc {X)p{x)dx « 2^ V loc.j Pjhdens 

■ d / 2 + RS j=-nN core 

Ncore 

= Yl w loc,J P J h i 

J— Ncore 



(9) 



where 



s \ ^ l^nj+k -A-J\ hard,s / 1f .\ 

W loc,J — 2-^1 ZJ, V loc,nJ+k' l iU ) 

k=—n 



v^t' s = v i^c d ' S ( x j ~ anc ^ P J * s ^ e electron density on the coarse-grid point Xj. 

As shown in Eqs.flZJ) and Q, the inner products have been replaced with the summation 
over coarse- grid points inside the core region, which produces only a modest overhead in the 
computational cost. Note that the weight factors w s j arising from the interpolation are inde- 
pendent of the wave functions, but dependent only on the known values of pseudopotentials 
on dense-grid points. Thus, if once computing the factors w s j for every molecular-dynamics 
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time step, we do not have to recalculate them throughout self-consistent iteration steps. 
The extension of the above procedure to the cases of higher-order interpolations is straight- 
forward. 



B. Illustration of Double-Grid Efficiency 

We now examine the efficiency of the timesaving double-grid technique by incorporating it 
into the formalism of the RSFD method. Hereafter, we obey the nine-point finite-difference 
formula, i.e., Nf = 4 in Eq. (jAlj) . for the differentiation of the wave function. The dense-grid 
spacing is set as h d ^ ns = h^/3, where (/i = x, y, and z) is the coarse-grid spacing in the \x 
direction. The ninth-order Lagrangian interpolation is used in the dense-grid interpolation. 
The exchange-correlation effects are treated as the local-spin-density approximation ^| 
within the framework of the density-functional theory j^D, 21]. The pseudopotential database 
NCPS97 [22! based on the pseudopotential of Troullier and Martins [3] is used for the ionic 
potential. The isolated boundary condition is imposed for all directions of x, y and z. 



wdrogen atom function of the cutoff 
I, we defined an equivalent energy cutoff 



First, the convergence of total energy for a 
energy is depicted in Fig. 3. According to Ref. 
E c ° ars [= (ir/h^) 2 Ry] to be equal to that of the plane-wave method which uses a fast-Fourier- 
transform grid with the same spacing as the present calculation. Hydrogen is one of the 
most difficult atoms to treat by the RSFD method because the s component of its nonlocal 
pseudopotential oscillates in the vicinity of the nucleus. The hydrogen atom is placed in 
the center of the neighboring grid points for the x, y and z directions. Without the use of 
the double-grid technique, the total energy does not converge even when the cutoff energy 
increases to 17 Ry. On the other hand, when the double-grid technique is employed, the 
total energy converges even with a small cutoff energy; the total energy converges rapidly 
and monotonically as the cutoff energy increases. 

Next, Fig. 4 shows the total-energy variation as a function of the displacement of the 
fluorine atom relative to coarse-grid points along a coordinate axis. The coarse-grid spacing 
is set at 0.33 a.u. The calculation of a model containing a halogen atom requires a quite high 
cutoff energy because the pseudopotentials of halogen atoms vary rapidly within the core 
region. Particularly, in the calculation using a discrete grid, such as the RSFD method, the 
total energy changes in an unphysical manner depending on the relative position of the grid 
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point and the nucleus when the grid spacing is coarse, which prevents us from implementing 
practical calculations. As seen in Fig. 4, the energy variation in our scheme is ~0.01% of 
the total energy (650.4 eV), which is negligibly small. 

As the third example, we apply the timesaving double-grid technique to the calculation of 
the equilibrium bond length of F 2 and Cu 2 molecules in order to examine the performance of 
the double-grid method for molecules. The coarse grid spacings are set at 0.33 a.u. for the F 2 
molecule, and 0.27 a.u. for the Cu 2 molecule. Figure 5 shows the adiabatic potential curves 
for the respective molecules. The results without any interpolation have many "humps", 
and the distances between adjacent humps on the respective dotted curves are close to 
twice the length of the grid spacing h^, which confirms that the oscillation depends on 



the relative position o 
problem, Gygi et al. 



t 



re atom with respect to coarse-grid points. To circumvent this 
| proposed a real-space grid in adaptive curvilinear coordinates. 
However, owing to the use of these distorted coordinates, their scheme needs Pulay forces in 
molecular-dynamics simulations; this is computationally very demanding. On the contrary, 
our approach requires only the Hellmann-Feynman forces, which significantly reduces the 
computational cost concerning the calculation of forces. As shown in Fig. 5, the locations 
of the minima in the result using the double-grid technique agree with the equilibrium bond 
length (2.68 a.u. for F 2 and 4.19 a.u. for Cu 2 ) obtained by experiments and other methods. 
These results clarify that the present method solves this problem and is very efficient and 
applicable for practical simulations. 



III. FAST SOLVER FOR POISSON EQUATION 

A. Fuzzy cell decomposition and multipole expansion technique 

The Hartree potential in the Kohn-Sham equation is commonly evaluated by solving the 
Poisson equation 

V 2 v H (r) = -4vrp(r). (11) 

In the case of isolated boundary conditions, one needs to determine the boundary values 
of the Hartree potential just outside of the calculation domain to set up the matrix for 
solving the Poisson equation [24 1. For small systems, the numerical summation over the 
grid is the most direct procedure. On the other hand, when the system becomes large, the 
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numerical summation requires so many operations as to result in a demanding scheme, since 
the computational cost for the direct summation is proportioned to the 5/3 power of the 
model size. Chelikowsky et al. |25j, |26J] proposed the procedure using a multipole expansion 
of the charge density around an arbitrary point r°. 

p(r') 



1=0 



T dr' 



p{r>) 



P z (cos 9')dr' 



Jp(r')dr' | y (r»-r°) 

I <y»0 I ' J ^ ^ I 



013 



3 (»> - r®)(r v - r°) - 5^ u \r 

I r*r* fr* I 5 



,0|2 



+ 



Here, p and z/ denote x, y and z. In addition, the functions Pz(cos#') (/ = 0, 1,2, 
the Legendre polynomials and cos 9' is expressed as 



cosfl' 



is the dipole moment, 



i>, I l< - rl)p{r')dr', 
and is similar to quadrupole moment |27| . 

1 



K-rl){r' v -rl)p{r')dr'. 



(12) 



are 



(13) 



(14) 



(15) 



By virtue of this method, one can make the computational cost proportioned to the model 
size. However, the accuracy of the solution largely depends on the choice of the position r° 
because the expansion is carried out only around one point. 

We now introduce the fuzzy cell decomposition and multipole expansion (FCD-MPE) 
method, which is free from this problem. A weighting function u s (r) for the multiple-center 
system centered at the s-th nucleus is prepared. This function is the defining function of 
so-called Voronoi polyhedra Q s which provides Wigner-Seitz cells and is set to satisfy 
the following equations. 

5>.(r) = l (16) 



S 



,1 € fi s 

".(r) = { (17) 
otherwise 

The charge density is divided into the charge distribution existing around the s-th nucleus 
using the weighting function u s (r). 

Ps {r) = p{r)uj s {r), (18) 



p(r) = J2ps(r) (19) 



By performing the multipole expansion for each p s (r) centered around the position of each 
nucleus R s , we have 



s \ 1 1 /u 1 1 

s 3( rfl -R^(r u -Ri)-S^\r-R s \ 2 



|r-H s | 5 



+ 



(20) 



where p and z/ are x, y and 2. In addition, and are 

P» = J(r',-R;)p s (r')dr', (21) 

and 

= / \K - RS M - K)Ps{r')dr'. (22) 

Note that if u> s (r) were in the manner of a step function at the boundary of Q s , many 
terms would be required for the expansion of Eq. (J2(Jj) . In order to carry out the expansion 
with as small a number of terms as possible, the behavior of u. q (r) near the boundary should 
be made as smooth as possible using the fuzzy cell technique |29(, in which the section of 
the boundary of Eq. ()17j) is fuzzy. When the fuzzy cell is employed, the multipole expansion 
up to the quadrupole is sufficient to obtain an accurate solution. 
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B. Efficiency of fuzzy cell decomposition and multipole expansion method 



As a numerical test of the the FCD-MPE method in accuracy, the Poisson equation is 
solved employing the boundary values obtained by the FCD-MPE method and those by 
the multipole expansion around one point in the calculation domain. Here, the charge 
distribution is assumed to be 

p{r) = J>^ x (— ) l exp[- 72li (r-a i ) 2 ], (23) 

i=l,2 

which imitates the charge distribution of a diatomic molecule, where a± = (2,0,0), 0,2 = 
(-2, 0, 0), 7i,i = 6, 7i, 2 = 4, 72,1 = 0.8, and 72,2 = 0.6. The cell size is set at 18.0 x 16.0 x 16.0 
a.u., a nine-point finite difference formula, i.e., Nf = 4 in Eq. (jAljl . is adopted for the second- 
order derivative, and the grid spacing is chosen to be 0.15 a.u. which is smaller than those 
of practical calculations so as to suppress the numerical error due to the finite-difference 
approximation. In the case of using Eq. (|12jl. the expansion is implemented around 

r o = Sr'p(r')dr' 

J p{r')dr> ' v ; 

Figure 6 shows the relative errors between the values of the Hartree potential for the cross 
section at z — 0.075 a.u. using the computed boundary values and the analytical boundary 
value, i.e., 

( x erf(^72-(r-0) 
Mr) = 2^ 7M x | r - a .| ' (25) 

8=1,2 ' l ' 

As is evident from these results, the accuracy of the calculation is improved markedly by 
the FCD-MPE method. 



IV. KOHN-SHAM HAMILTONIAN UNDER TWIST BOUNDARY CONDITION 

A. Discretization of screw operator for twist boundary condition 

We here address the discretization of the Kohn-Sham equation under the twist bound- 
ary condition (see Fig. 7 for an example). The boundary- condition operators in x and y 
directions, t x and T y (refer to Appendix), are zero and the operator in z direction, T z , is 
represented by using a screw operator S j33j, which is an N xy {= N x x A y )-dimensional ma- 
trix in terms of a translation L z down (up) the z axis in conjunction with a right-handed 
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rotation tp (—(p) about the z axis. Here, N x and N y are the numbers of grid points in the x 
and y directions, respectively. 

Let us consider the discretization of the screw operator S. The wave function multiplied 
by the screw operator is defined as 

\V{z k )) = S{<p)\*{z k )) 

» S d (tp)\V(z k )), (26) 

where ^(z k ) is a columnar vector consisting of N xy values of the wave function on the x-y 
plane at the z = z k point, ip Zk (xi,yj). Note that the wave function ip' cannot be evaluated 
exactly in the real-space formalism with the Euclid coordinate, since the twist angle ip is 
not always in a multiple number of tt/2. Therefore, the wave function is obtained by an 
interpolation using the values of the wave function on the neighboring grid points. In the 
case of the linear interpolation, we have 

+^z k (xi>,yj'+i) + Uipz k {xv+i, yj'+i), 

(27) 

where 

6 = (xi'+i - Xi cos p> + yj sinp)/ h x 

x (y j/+ i - Xi sin tp - yj cos tp)/h y , 
^ 2 = (xi cos tp-yj sin tp- Xi>)/h x 

x (y j/+1 - Xi sin tp - yj cos <p)/h y , 
£ 3 = (x v+x - Xi cos p + yj sin tp)/h x 

x (xi sin tp + yj cos tp — yy)j h y , 
U = (xi cos p - yj sin tp-x v )/h x 

x(xiSin tp + yj cos tp — yj/)/h y . (28) 

Here, Xi>, av+i, yj> and yj'+i are so determined as to satisfy xy < Xi cos <p—yj sin tp < xy + i and 
Vj' < Xjsin<^ + yjcostp < yj>+i- Whereas the exact screw operator holds S(— tp) = S^(tp), 
the discretized screw operator does not necessarily satisfy the identity due to the loss of 
accuracy of the interpolation using discrete grid points. The boundary condition T z in 
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Eq. (jSl is defined as e ik * Lz [S d {<^) + «Sj(-^)]/2 so as to keep the discretized Hamiltonian as 
an Hermitian, where k z is the Bloch wave number. Concerning the long-range ionic potential 
in the case of the twist boundary condition, refer to Ref. j^. 



B. Application: Band structure of helical carbon nanotube 



Carbon nanotubes [30j consist of a few concentric tubes each of which has carbon-atom 
hexagons arranged in a helical fashion about the axis. Their topologies are characterized 
by a set of two numbers (n, m) called chiral indices. Since they were first discovered, their 
atomic and electronic structures have been intensively explored by both experimental and 
theoretical approaches. Previous theoretical studies have shown that the carbon nanotubes 
change between having metal and semiconductor properties depending on their structure 



such as their diameter and helical arrangement [3JJ, |32j. In order to present the efficiency 
and applicability of the RSFD method, we compute the electronic band structure of the 
(5, 5) carbon nanotube using the twist boundary condition. 

The computational conditions are as follows: the influences of core electrons are induced 



by the Troullier-Martins-type pseudopotential 



form 







22j in the Kleinman-Bylander nonlocal 



17| . The tube structure is generated with a radius of 6.41 a.u., a screw operation 
having a twist angle of tt/5 rad and a translational shift of L z =2.32 a.u. chosen to yield 
nearest-neighbor separations between rings equal to the in-ring values. Exchange-correlation 



effects are treated using the local density approximation [19| and the central finite-difference 
formula, i.e., Nf = 1 in Eq. (|A1|) . is adopted. The coarse grid spacing is set at 0.33 a.u. and 
the denser grid spacing is set at 0.11 a.u. in the vicinity of the nuclei with the augmentation 
of double-grid points. 

We show in Fig. 8(a) the helical band structure of the (5, 5) carbon nanotube. Since there 
are multiple equivalent ways of depicting the band structure of (n, n) carbon nanotubes, 



Mintmire et al. 



341 ] proposed a pseudoband structure. By their procedure, we can classify 



all the bands into two representations; an "a" representation which varies in the Brillouin 
zone depending on possible helical operators and an "e" representation which moves in the 
Brillouin zone depending on helical twist angle. The pseudoband structure is obtained by 
shifting all the "e" bands in k by a quantity corresponding to a phase factor of ±j(p. This 
choice of phase shifts the "e" bands to doubly degenerate bands. 
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The pseudoband structure of the (5, 5) carbon naotube is illustrated in Fig. 8(b). We also 
depict in Fig. 8(c) the band structure using the conventional supercell which has twice the 
length in the tube direction. One can recognize that the present pseudoband structure leads 
to a depiction of the bands as equivalent as possible to the bands that would be present if 
the system had translational periodicity with the repeat distance of the helical supercell. 

V. SUMMARY 

We have presented efficient techniques for the RSFD calculations together with several 
examples. The main points of our method are as follows: (i) The timesaving double-grid 
method can suppress the spurious oscillation of the total energy for the displacement of 
the atom relative to the coarse-grid points, and there is little computational effort for the 
double-grid method by virtue of the integration over the coarse-grid points, (ii) The FCD- 
MPE method improves the accuracy of the computation for the boundary values of the 
Hartree potential required in the case of the isolated boundary condition, (iii) We can easily 
treat a system having a twist angle at its boundary by including a twist operator in the 
Kohn-Sham Hamiltonian. The pseudoband structure of the (5, 5) carbon nanotube agrees 
with that calculated employing the large supercell. 

From what has been discussed above, our techniques enable us to implement highly accu- 
rate first-principles calculations based on the RSFD approach and are efficiently applicable 
to various systems, in particular, a system having a twist angle at its boundary. When 
coupled with conductance calculations, this approach would be useful in studying electron 
transport behavior in helical nanotubes P, [J. 
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APPENDIX A: KOHN-SHAM EQUATION IN REAL-SPACE FINITE- 
DIFFERENCE FORMALISM 



We demonstrate the procedure for discretizing the Kohn-Sham Hamiltonian 20|, |2l| in 
the real-space formalism. In the conventional methods using basis sets, the derivative of 
wave functions arising from the kinetic-energy operator in the Kohn-Sham Hamiltonian can 
be computed by differentiating the basis. However, since the prime concept of the RSFD 
method does not employ any bases, the second derivative of function f(x) at grid point 
x = ih x (i: integer) is approximated by the following finite-difference formula. 



d 2 



N f 



x=ih-r 



^ Cnf(ih x + nh x ) 



(Al) 



n=-Nt 



Here, Nf and h x are the parameters determining the order of the finite-difference approx- 
imation and the grid spacing, respectively. In the case of Nf = 1, the coefficients are 
c_i = ci = l/h x , and c = —2/h 2 x . See Ref. j3] for the values of coefficients c n when 
Nf > 1. Although we here present the case of the central finite difference, i.e., Nf = 1, 
and local pseudopotentials in order to describe the essence of this scheme, the inclusions of 
higher-order finite-difference formulas and nonlocal parts of the norm-conserving pseudopo- 
tentials are straightforward. 

In the RSFD method, the wave functions and electron density are given only on discrete 



grid points in three-dimensional real space. Therefore, the Kohn-Sham Hamiltonian [21 1 
acting on a wave function must be given in a form discretized in real space. The discretized 
Kohn-Sham equation is written as 



V{z x ) B z 
B\ V(z 2 ) B z 
'•• '•• '•• 



B\ V(z k ) B z 











B\ V(z Nz -i) 
B\ 



z 







B z 
V{z Nz ) 





tt(*l) 


















= E 





















(A2) 



where T z is an N xy (= N x x iV y )-dimensional matrix specified by the boundary condition 
in the z direction [for example, see Eq. (|A9)) ]. ^(zk) is a columnar vector consisting of N xy 
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values of the wave function on the x-y plane at the z = Zk point, ip Zk (xi, i/j), B z is a constant 
matrix proportional to TV^-dimensional unit matrix J, 

1 



Bz= 2hl 



I, 



(A3) 



and N x , N y and N z are the numbers of grid points in the x, y and z directions, respectively. 
In addition, V(zk) is the following N xy - dimensional matrix defined on the x-y plane at the 
z = Zk point. 



V{z k ) 



V Zk ( Vl ) B y ... 

£t V Zk (y 2 ) B y 

'•• '•• '•• '•• 

'■■ '•• Bl V Zk (yj) By 







v 



T 







B\ V Zk (y Ny ^) B y 
Bl V Zk (y Ny ) 



where B y is an A^-dimensional matrix similar to B z , 



y 2hl 



I. 



(A4) 



(A5) 



and the block matrix V Zk (yj) is defined on the x line at the (yj, z k ) point as 



v yj , Zk (xi) b x 



X 





v*(vs) = 



v yj ,z k (x 2 ) b x 



X 





(A6) 



with 



and 



K v Vi^k( X N x -l) 

6! 



v yi , Zk (xN x ) 



h 1 1 1 



, N 1 1 1 \ 

v yj , Zk {Xi) = j^ + j^ + j^ + v{x i: y j: z k ). 



(A7) 
(A8) 
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Here, T y (t x ) is an A^-dimensional matrix (a coefficient) determining the boundary condition 
in the y (x) direction and v is the sum of the external, Hartree, and exchange-correlation po- 
tentials. In the case when the periodic boundary conditions are imposed on all the directions, 

T z = e ikzLz B z , (A9) 

T y = e ikyLy B y , (A10) 

and 

r x = e ik * L °b x , (All) 

where k x , k y , and k z are the Bloch wave numbers, and L x , L y , and L z are the unit-cell 
lengths. On the other hand, when the isolated boundary conditions are imposed, T z = 0, 
T y = 0, and r x = 0. 
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FIGURE CAPTIONS 



FIG 1: Double-grid adopted in the text. The "x" and "•" correspond to coarse- and 
dense-grid points, respectively. The circle shows the core region of an atom which is taken 
to be sufficiently large to contain the cutoff region of nonlocal pseudopotentials. 

FIG 2: Functions on coarse- and dense-grid points in the one-dimensional case. Xj (xj) 
represents a coarse(dense)-grid point with j = nJ + 1 (0 < t < n), and so Xj = x n j. (a) 
The inner product between the wave function i]){x) and the nonlocal part of pseudopotential 
v(x) evaluated over coarse-grid points, (b) The inner product evaluated over dense-grid 
points. 

FIG 3: Convergence of the total energy for the hydrogen atom as a function of the 
coarse-grid cutoff energy E c ° ars . The atom is located at the center between adjacent 
coarse-grid points. 

FIG 4: Variation of the total energy as a function of the displacement of the flourine atom 
relative to coarse-grid points along a coordinate axis. The coarse-grid spacing is 0.33 
a.u. 

FIG 5: Relative positions with respect to the grids for F 2 and Cu 2 molecules and adiabatic 
potential curves, (a) The center of gravity of the molecule is placed at the center of 
neighboring grids with respect to the x, y and z directions. In the figure, " x " indicates the 
coarse grid and • represents nuclei, (b) The center of gravity of the molecule is placed on 
the grid plane which orthogonally intersects the coupling axis. For the grid plane parallel 
to the coupling axis, the center of gravity is placed at the center of neighboring grids, (c) 
Adiabatic potential curves of F 2 molecules. The grid spacing is 0.33 a.u. (d) Adiabatic 
potential curve for Cu 2 molecules. The grid spacing is 0.27 a.u. The calculated points are 
fit to spline functions as a guide to the eye. 

FIG 6: Relative errors between the Hartree potential obtained using boundary values 
computed by each multipole expansion method and that evaluated by the analytical 
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solution for the cross section at z =0.075 a.u. Results using boundary values (a) obtained 
by the FCD-MPE method and (b) obtained by the multipole expansion around one point. 
Spheres show the locations of the nuclei. 

FIG 7: (a) Front and (b) axial views of the (5, 5) carbon nanotube. r, ip, L z are the radius, 
the helical twist angle and the translational shift of the nanotube, respectively. 

FIG 8: (a) Electronic band structure of the (5, 5) carbon nanotube using the supercell with 
length L z , (b) pseudoband structure, and (c) electronic band structure using the conventional 
supercell with the length 2L Z . The zero of energy is chosen to be the Fermi level. 
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